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Experiments on particles' motion in living cells show that it is often subdiffusive. This subdiffusion 
may be due to trapping, percolation-like structures, or viscoelatic behavior of the medium. While the 
models based on trapping (leading to continuous-time random walks) can easily be distinguished 
from the rest by testing their non-ergodicity, the latter two cases are harder to distinguish. We 
propose a statistical test for distinguishing between these two based on the space-filling properties 
of trajectories, and prove its feasibility and specificity using synthetic data. We moreover present a 
flow-chart for making a decision on a type of subdiffusion for a broader class of models. 



Experiments on particles' motion in living cells aimed 
on understanding molecular crowding [1-4] have unveiled 
that diffusion in such environments is often anomalous, 
i.e. the mean squared displacement (MSD) does not grow 
proportionally to time, (x 2 (t)) oc t , but follows a slower 
pattern 

(x 2 (t))^e (i) 

with < a < 1 (subdiffusion), and the nature of this 
anomaly has to be understood. Anomalous diffusion is 
not only apparent in biological systems, but is found in 
complex systems ranging from amorphous semiconduc- 
tors [5] , goeology [6] , to turbulent systems [7] . 

There are several mathematical models leading to sub- 
diffusion, corresponding to different physical assumptions 
about the structure and energy landscape of the system 
in which the subdiffusive motion takes place. Since one 
is mostly interested in the actual microscopic structure 
of the system, an important task is working out tests 
which allow for distinguishing between different models 
giving the same prediction for the MSD. The three most 
popular models which might be pertinent to explaining 
subdiffusion in cells are: 

(i) continuous time random walk (CTRW), a math- 
ematical model which is physically realized in systems 
with traps, i.e. binding sites of different energetic depths, 
a case pertinent to energetic disorder, 

(ii) diffusion on fractal structures, as exemplified by 
percolation, a situation pertinent to structural disorder, 
and 

(iii) fractional Brownian motion [8] (fBm), a Gaussian 
process with stationary increments which satisfies the fol- 
lowing statistical properties: the process is symmetric, 
i.e {xnifii) = where xh(0) = 0, and the MSD scales as 
{x 2 H (t)) ~ t 2H where H is the Hurst exponent. Note that 
H < 1/2 leads to subdiffusion, while H = 1/2 recovers 
Brownian motion. fBm physically corresponding to sys- 
tems with predominating slow modes of motion and is 
realized in viscoelastic media as exemplified by polymers 
and polymer networks, where disorder does not play a 
leading role. 



Lastly, one has to discuss 

(iv) the time-dependent diffusion coefficient (TDDC) 
model - normal diffusion with a time-dependent diffu- 
sion coefficient, which is used to fit experimental re- 
sults from, for example, FRAP (fluorescence recovery af- 
ter photobleaching) experiments [9]. This model corre- 
sponds e.g. to a situation when the step rate is explicitly 
time-dependent, and does not have a clear physical in- 
terpretation in application to crowded media. 

The non-stationary (and non-ergodic) models of 
anomalous diffusion like CTRW or TDDC are easily dis- 
tinguished from the ergodic and stationary models of dif- 
fusion (as exemplified by fBm or diffusion on percolation 
structures) by applying tests aimed onto checking sta- 
tionarity of increments or ergodicity. At present, two 
of them can be recommended: the p-variation test [10] 
which can be considered as a generalized test of tempo- 
ral homogeneity of the process, and the moving average 
vs. ensemble average test [11] which is a clear test for 
ergodicity, see [12, 13] for their practical application. 

It is much harder to distinguish within the class of er- 
godic non-Markovian processes, i.e. to tell whether the 
observed subdiffusion is due to geometrical restrictions 
(e.g. percolation) or to a viscoelastic medium (fBm). 
More detailed information on these two models, and how 
to simulate them, appears in the Supplementary Mate- 
rial. Note that both models correspond to antipersis- 
tent random walks (RWs), and may have the same step- 
step (or velocity-velocity) correlation function. The cor- 
responding correlation function for percolation is calcu- 
lated in [14] and on the coarse-grained level it is con- 
nected with the spectral dimension of the percolation 
structure. The step-step correlation functions are shown 
in the Supplementary Material. For any physical model 
resulting in fBm the correlation function follows from the 
spectral properties of slowest modes. Distinguishing be- 
tween the models is particularly challenging in single- 
molecule tracking experiments where only one or few tra- 
jectories of motion are recorded [15]. 

One fundamental difference between the two is the 
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probability distribution function (PDF) of displacements 
which is Gaussian for fBm but non-Gaussian for perco- 
lation, meaning that a Gaussianity test (i.e. in the exact 
relation between the second and the higher even central 
moments of the PDF) may in principle solve the prob- 
lem [16]. However, the limited amount of available infor- 
mation is not enough to produce a distinguishable PDF 
(an example is shown in the Supplementary Material), 
and moreover no analytical form is known for the perco- 
lation PDF. Moreover, Gaussianity on its own does not 
shed light on the nature of the type of constraint govern- 
ing the tracer's motion, i.e whether its motion is confined 
to an inhomogcneous geometrical structure, which does 
not considerably change on the time scale of the experi- 
ment, or such a structure is absent, and the restrictions 
to the motion change with time (like in the Rouse model 
of polymers or in single-file diffusion) . This information 
can be delivered by the tests of spatial homogeneity of 
the corresponding motion. The aim of the present work 
is to give such a test on a single trajectory level, and to 
prove its feasibility and specificity using synthetic data 
for percolation and for fBm with exactly the same MSD 
behavior. 

Our present discussion is confined to a two-dimensional 
(2d) situation, such as the diffusion of membrane proteins 
in the cell membrane, which constitutes one of the most 
interesting cases where single molecule tracking methods 
are used, see e.g. [17] and references therein. Our dis- 
cussion can easily be generalized to 3d, if the data for 
all three coordinates are available. Caution is advised if 
the data available corresponds to the 2d projection of a 
3d trajectory, like in [18], in which case our method may 
not be appropriate. 

On the level of the RW description, the processes with 
non-stationary and with stationary increments differ in 
how the clock time t is translated into the steps of the 
problem. In both CTRW and TDDC the steps follow 
inhomogeneously in time, and the mean number of steps 
n taken up to time t grows as (n) oc t a , while the MSD 
as a function of the numer of steps grows as (x 2 (t)} oc n. 
Thus, CTRW and TDDC models correspond to normal 
diffusion if the clock time is translated into steps of the 
RW process. This transformation can cither follow a 
random process (CTRW) or be deterministic (TDDC). 
These processes fill space homogenously like in normal 
diffusion. On the other hand, for fBm and for a RW on 
a percolation cluster, being time-homogeneous processes 
(with stationary increments), the number of steps is al- 
ways proportional to time. Here the fractal dimension of 
the trajectory is connected to the exponent a character- 
izing anomalous diffusion, 

(r 2 (t)} oct" = t 2 ' d ™ «n 2 ^, (2) 

where d w = 2/ a is called the walk dimension. For fBm 
the exponent a is related to the Hurst exponent by a = 
2/d w = 2_ff.The fractal dimension df is defined through 



how the amount of available sites within a radius r scales 
with r: M n ~ r d f . 

Let us consider the number of sites within a radius r 
as a function of the average time needed to reach such 
a radius. We do so by substituting the square root of 
the MSD's time dependance in place of r: M n ~ r df <~ 
(„!/*.)*/ = n <V2. 

Note that we are dealing with a recurrent walk, where 
M n grows slower than n, i.e df < d w .ln this case each 
of the sites within the reachable distance is visited at 
least once, and the total number of distinct visited sites 
behaves as S n w M n , i.e [21, 22]: 

S n ~n d i' d ™. (3) 

In the case of fBm the geometry is Euclidean, meaning 
that df = 2 or 3 in 2d or 3d respectively. In percolation, 
on the other hand, df s=s 1.8958 and df w 2.52 when 
embedded in 2d and 3d respectively. The walk dimension 
associated with a RW on a 2d percolation cluster is d w — 
2.87. 

In our approach to the problem we propose to exploit 
the fact that fBm explores an Euclidean structure (df — 
d, with d being the dimension of space), while a RW on 
a percolation cluster explores a fractal one with df < d. 

To formulate the null-hypothesis as to how the RWcr 
fills space, let us look at the ratio of the average number 
of distinct visited sites within n time steps, S n , and the 
space enclosed in the radius r(n) which the RWer reaches 
on average within the same number of time steps, r d (n) = 

(r 2 (n)) d / 2 . We examine S n /(r 2 (n)) d / 2 ~ n^r , so that 
our test is based on the calculation of the exponent: 



i.e. the difference between df and d for a given d w . We 
note that if the RWer fills space homogeneously, the two 
quantities grow at the same rate, meaning that the curve 
is expected to be flat, or 5 = 0. Indeed for fBm df = d, 
meaning that 5 — 0, as opposed to the case of a RW on 
a percolation cluster where df < d, leading to 5 < 0. 

To assess the success of our test, we simulated 2d single 
trajectories of fBm and of a RW on a percolation cluster 
at criticality. We chose the Hurst exponent H so that 
the MSD of the two is identical; the value d w = 2.87 for 
a RW on a percolation cluster corresponds to H = 0.348. 
We modeled an experiment with optical limitations, us- 
ing a coarse grained lattice with a characteristic grain 
size A. Two sample trajectories are shown in the inset of 
Fig. 1. Note that the trajectory of the RW on a per- 
colation cluster is restricted here to only horizontal and 
vertical directions since the percolation cluster is based 
on a square lattice (see Supplementary Material for an 
example of such a cluter). This is also the reason for the 
small oscillations found in the ACF, also shown in the 
Supplementary Material. A fast and precise generator 
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FIG. 1: (color online) Participation function divided by the 
MSD (S n /(r 2 (n)) d/2 for 2d), temporally averaged with a mov- 
ing window of < r < 150 for 5 trajectories of fBm (blue, 
which flatten out) and 5 trajectories created of a RW on a 
percolation cluster (red, with a clear negative slope). All tra- 
jectories are 40000 time steps long. Two sample trajectories 
prior to coarse graining are shown in the inset. 

for fractional Gaussian noise in the antipersistent case is 
described in [23]. 

Thus, our test is based on calculating 8 from the slope 
of S(t)/(r 2 (t)) d / 2 on the double logarithmic scale and 
testing whether this S is different from zero. The null- 
hypothesis 5 = corresponds to fBm, and its rejection 
witnesses in favor of the percolation model. For the spe- 
cific case of a RW on a 2d percolation cluster, we expect 
8 = —0.037. Note that it is typically hard to detect the 
differences in the exponents of such magnitude on the 
basis of relatively short runs. However, as will be seen in 
what follows, we are in luck. 

We found that for single trajectories as is the method is 
not sensitive enough due to strong noise. This noise can 
be reduced by looking at a moving-window time average. 
Fig. 1 displays temporally averaged S{t) / {r 2 {t)) d / 2 in 2d, 
for trajectories of a RW on a percolation cluster, and of 
fBm. The fBm curves flatten after the first couple of 
steps, as expected, whilst the percolation ones clearly 
have a negative slope. 

We now fit each of these curves to a power-law and ex- 
tract the exponent corresponding to 5. Fig. 2 shows the 
distribution of i5 resulting from 400 fBm and percolation 
trajectories. A closer look at Fig. 2 reveals that whilst the 
peak of the fBm 8 distribution is centered around 0, the 
percolation distribution is not centered around —0.037, 
but at a much larger negative number: i.e. around —0.18 
for an averaging time window T max = 50 and around 
—0.12 for T max = 550. This is due to large corrections 
to scaling for the percolation case (see e.g. [24]), which 
luckily play in our favour: for smaller T max the distribu- 
tions are clearly distinguishable, with no overlap, mean- 
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FIG. 2: (color online) Distribution of <5 for 400 fBm trajecto- 
ries (blue peaks on the right) and 400 trajectories of RW on 
a percolation cluster (red peaks on the left). The trajectories 
are temporally averaged within time windows of T max = 50 
time steps (solid line) and of T max = 550 (dashed line). While 
the fBm distributions stay centered at the expected value 
of 0, the ones for percolation start far off (at w —0.18 for 
T m ax = 50) and slowly converge to the expected value of 
—0.037 as shown in the inset. The dashed line indicates the 
expected asymptotic value. 



ing also that a relatively short trajectory is enough. So 
in practice, given a single trajectory one may calculate 
8 for different T max , and see whether these are negative 
and converge to a smaller negative number (pointing at 
a RW on a fractal), or tend to zero (fBm). We can now 
add this test to the test of ergodicity, building a toolbox 
to help identify the underlying physics of a given experi- 
mental trajectory. We summarize our toolbox in the form 
of a decision tree in Fig. 3. 

One may take this toolbox one step further and con- 
sider a more general and realistic scenario, of subordi- 
nation of any two of these models, as previously consid- 
ered [25, 26]. In a biological cell, for example, there is 
no reason why one may not encounter both energy traps 
(modeled with CTRW) and crowding (modeled as a RW 
on a percolation cluster), i.e the problem would be mod- 
eled as a RW on a percolation cluster, with the subordi- 
nation of waiting times at each step. This generalization 
is out of the scope of this paper, and will be set forth 
elsewhere. 

We proposed here a toolbox of tests that may be run on 
single trajectories in the aim of discerning between pos- 
sible physical realities including combinations of energy 
traps, structural disorder or crowding, and a viscoelastic 
medium. We note that not all tests may be feasible ac- 
cording to the experimental setup and the type of data at 
hand, but it nontheless illuminates the possibilities and 
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FIG. 3: (color online) Flow chart of the decision process in dis- 
cerning between the three main subdiffusive models: CTRW, 
fBm and a RW on a fractal structure. One starts by assessing 
the ergodicity or temporal homogeneity of the process. If the 
process if found to be non-ergodic, it is CTRW. If the process 
is ergodic, one is left to discriminate between fBm and a RW 
on a fractal structure. Here one analyses S(t)/(r 2 (t)) d ^ 2 for 
2d. If this ratio is a constant, the process is fBm, if it decays, 
the process is a RW on a fractal structure. 



gives a broader understanding. 
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DESCRIPTION OF THE MODELS 



delbrot and van Ness [3] denned fBm as: 
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Let us first describe how one simulates a percolaton 
cluster at criticalitiy [1]. We take a two dimensional lat- 
tice, and start removing sites. The critical concentration 
of sites for a two-dimensional square lattice is p c = 0.5. 
For a higher concentration there are multiple ways for a 
RW to reach from one side of the lattice to the other, for 
a lower concentration so many sites have been removed 
that the two sides of the lattice are detached. The critical 
concentration is that in which, on average, removing one 
more site will detach the lattice sides. A RWer is allowed 
to walk between the remaining sites, with the remaining 
bonds. Fig. 1 shows and example of such a percolation 
cluster at criticality. Note the empty spaces at differ- 
ent (scale-free) sizes. It should also be noted that since 
the cluster is based on a square lattice, the trajectory of 
a RW on such a cluster will also be restricted to only 
horizontal and vertical steps. 

FBm is generated from fractional Gaussian noise, just 
as Brownian motion is generated from white noise. Man- 
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where T is the Gamma functin and < H < 1 is the 
Hurst parameter. Note that H = | recovers ordinary 
Brownian motion. 

As for fBm, quite a few methods of simulation exist. 
We chose one of the most popular methods which uses 
a discrete Fourier transform. The method is thoroughly 
described in [2]. 



PDF AND ACF ARE INADEQUATE TO 
DISTINGUISH BETWEEN THE MODELS 



Fig. 2 and 3 show the PDFs for two pairs of trajectories 
of fBm and for RW on a percolation cluster. The behavior 
of these two examples is drastically different, and shows 
that obtaining the known final forms of the distribution 
needs large ensembles. Fig. 4 shows the ACF for fBm 
and for a RW on a percolation cluster. The two are iden- 
tical apart from very small oscillations apparent in the 
RW on a percolation cluster. Adding exponentially dis- 
tributed waiting times to the steps (a more realistic case 
than equally timed steps), eliminates these oscillations, 
yielding undistinguishable ACFs. 
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FIG. 1: An example of a two-dimensional site percolcation 
cluster at criticality. 
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FIG. 2: (color online) Probability distribution function calcu- 
lated form a single trajectory of fBm (blue) and of a RW on 
a percolation cluster (red). Both trajectories are 40000 time 
steps long. 
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FIG. 3: (color online) (color online) Probability distribution 
function calculated form a single trajectory of fBm (blue) and 
of a RW on a percolation cluster (red). Both trajectories are 
40000 time steps long. 
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FIG. 4: (color online) The step-step ACF of fBm (blue) and 
of a RW on a percolation cluster (red). In the inset: the ACF 
of the two but after adding waiting times from an exponen- 
tial distribution (a more realistic view than identically timed 
steps). The slight oscillations of the RW on a percolation 
cluster disappear, and the two ACFs are indistinguishable. 



